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1. OVERVIEW 

This grant concerns the development of embedded function methods to accurately and 
efficiently predict heat transfer and skin friction coefficients in high-speed compressible 
turbulent flows near walls, without computing the flow all the way to the wall. In general, the 
methodology being developed can be applied for a wide variety of outer region turbulence 
models, including two equation models. However, in the current program the effectiveness of 
the technique is being demonstrated using simple outer-region models similar to the Baldwin- 
Lomax and Cebeci-Smith eddy viscosity models. It is worthwhile to note that the embedded 
function method generally applies to attached turbulent flows where a turbulent wall structure 
exists and is attached to the wall. Such situations normally occur over a large portion of a 
computational domain. Virtually all outer region turbulent models are tied tightly to a 
logarithmic behavior in the mean profile near the wall. Experimental data for turbulent 
separated flows shows that the logarithmic behavior in the mean profile disappears as 
separation is approached; although some experiments indicate that a logarithmic profile can be 
established downstream in a region of backflow, the “constants” in the logarithmic zone are 
radically different from those in the upstream boundary layer. This is simply an indication 
that the physics near the wall in a separated turbulent flow are substantially different from 
those in an attached turbulent flow. For this reason, the present program has focussed on 
developing methods and testing of the concepts for attached flows rather than trying to force a 
conventional “law of the wall” into a zone of backflow; not surprisingly, previous attempts at 
the latter approach (using two-equation models) have not met with any significant degree of 
success (Rubesin and Viegas, 1985). Essentially, this is because although the dynamics of the 
near-wall flow in an attached turbulent boundary layer are relatively well documented (Walker 
et al, 1989), the dynamical features of a zone of reversed turbulent flow are not well 
understood and have not been thoroughly investigated experimentally. It is likely that an 
appropriate modification of the embedded function approach can be made for reversed flow at 
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a later stage. However, because the compressibility introduces effects and issues that have 
been dealt with only marginally in the literature, the present work has principally focussed on 
attached high-speed flows. 

The research program started in January, 1988, as a cooperative effort involving 
Lehigh University and United Technologies Research Center. The theoretical analysis and 
algorithms are being developed at Lehigh University under the direction of Professor J. D. A. 
Walker; the work forms the basis of the dissertation of Ph.D. candidate Mr, Jun He, who is 
supported as a research assistant on the grant. Implementation of the algorithms was to be 
carried out at UTRC by Mr. Greg Power in the NASA Langley code NASCRIN; this is an 
explicit Navier-Stokes code which utilizes the Baldwinn- Lomax turbulence model and uses a 
mixing length formulation (with a Van Driest damping factor) in order to calculate the flow all 
the way to the wall. Due to heavy contractual commitments at UTRC during 1988, Mr. 
Power was not able to devote the necessary time to effectively implement the embedded 
function methodology in the NASCRIN code; due to budget limitations, NASA support for the 
UTRC portion of the effort was not renewed for the second year (1989) of the project. In 
March, 1989, UTRC hired Dr. Seyhan Ersoy to carry out introduction of the embedded 
functions into NASCRIN; the unused funds from the first year of the grant as well as 
significant additional financial support from UTRC is being used to pay Dr. Ersoy’s salary 
(until March, 1990). 

At present the theoretical work at Lehigh is presently on the course laid out in the 
original research proposal for a three year effort. The wall function method has been extended 
up through the supersonic range to hypersonic speeds. The algorithms have been successfully 
introduced into the NASCRIN code at UTRC; testing of the procedures is currently being 
carried out for progressively more complex flow situations. 

2. OVERVIEW OF LEHIGH THEORETICAL WORK 

There are a number of fundamental issues relating to compressible turbulent flow 
which have not been adequately addressed in the previous literature. Unfortunately, there has 
been a trend to extrapolate turbulence models and modeling concepts for incompressible flows 
(where data is abundant) to the high-speed compressible case, essentially without modification. 
In order to ensure good success with an embedded function approach, it is vital to isolate the 
correct scaling laws in the overlap region. A number of the fundamental issues have been 
resolved by the present research. 
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A first question concerns the precise form of the law of the wall for compressible flow. 
An extrapolation of the incompressible law of the wall would suggest 



(i) 


in the overlap zone 1 . Here k is the Von Karman constant and = 5; in addition A t is a 
parameter proportional to the wall-layer thickness whose value is (y/u T ) for incompressible 
flow but whose form for compressible flow needs to be determined. The present research has 
shown that equation (1) is not adequate for compressible flows and that the “compressible law 
of the wall” is of the form 



( 2 ) 


where Y is the Dorodnitsyn variable 



( 3 ) 


Thus the normal variable in equation (2) implicitly contains the density; a form of 
A # - = jj w / u r Re is indicated by the analysis (Re here is the Reynolds number and the subscript 
w indicates a quantity is evaluated at the wall). 

A second important issue relates to the temperature distribution, since an important 
objective of the present analysis is to develop the embedded function method for flows with 
heat transfer. A central question is whether it is the static temperature or the total enthalpy 
which is logarithmic in the overlap zone. Unfortunately, there is very little data available for 
direct measurements of either static temperature or total enthalpy, particularly near the 
surface. The present analysis, as well as some very recent data confirms that the total 

1 Note that all variables in this report are dimensionless with the exception of those 
denoted with an asterisk superscript. 
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enthalpy is logarithmic and in the Dorodnitsyn variable; the general form is 



( 4 ) 


Here I„, is the total temperature at the wall, Pr is the Prandtl number and is the analogous 
quantity to the von Karman constant; in addition Q is a dimensionless heat transfer parameter 
at the wall defined by 


Q = 


-kt(OT7«y*) y*=o 

p*w u;h: 


(5) 


Here the subscript c denotes a quantity evaluated at the mainstream. 

A third issue relates to the common belief that incompressible turbulence models can 
be used without modification as the Mach number increases. The current research has shown 
that this cannot be done without incurring substantial error. For example, the outer region 
Baldwin-Lomax model is an eddy viscosity formulation with 


-pu'v' = pe( y) 


( 6 ) 


for which 


«(y) = 


t 


CpKY maiFmai 
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y > y m, 
y < y m. 


(7) 


where F ma x is the maximum of the function y|<9u/<9y| and y m ax is the corresponding ordinate; 
y m denotes the point where both portions of the function are the same. The model is modified 
in the wall layer to become a mixing length model with a Von Driest damping factor. The 
present research has shown that models like equation (6) are not self-consistent in a high-speed 
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compressible flow and will produce progressively worse results as the Mach number increases. 

Consider, for example, the eddy viscosity function. The wall layer is a region of 
constant total stress to leading order and the asymptotic analysis shows that 


-pu'v' — ► P.U? as Y + — ► oo. 


( 8 ) 


Here p w is the wall density and Y + is the scaled Howarth-Dorodnitsyn variable 


y+ = XjzM*, Y = £ pdy. (9) 

If the turbulence model used in the outer region is an eddy viscosity model, it will have the 
following form 


* p u'v' 



( 10 ) 


where e is the eddy viscosity function. In the outer layer, u is in a defect form which behaves 
logarithmically in the overlap zone (in the Howarth-Dorodnitsyn variable), viz. 


u ~ U e + u r 





(ii) 


as Y/A o— ► 0. Combining equations (9), (10), and (11), it follows that 
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( 12 ) 


This is the critical relation that a self-consistent turbulence model should satisfy in a 
compressible flow. With increasing distance from the wall, the eddy viscosity function must 
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deviate from the linear behavior in the overlap zone, and for low Mach numbers typical far- 
field eddy viscosity formulae are of the form 


e - K\Jc&\ 


(13) 


for the Cebeci-Smith model and 
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F 
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(14) 


for the Baldwin-Lomax model. During the past year, a number of models were considered 
with the requirements that: (1) c exhibit the asymptotic form (12) near the overlap zone; (2) 
for low M e , the model should reduce to either of equations (13) or (14); and (3) no additional 
correlations should be introduced. Unfortunately, there is a large number of possible models 
that could be considered which satisfy these criteria and several models were considered in the 
following way. 

A considerable data base (at least for the velocity profiles) exists where the flow is 
close to being in equilibrium. In order to test the embedded function procedures and a variety 
of potential models, the limiting case of self-similar flow was considered; in this limit the 
governing equations reduce to ordinary differential equations which can be solved accurately 
and quickly to facilitate comparison with measured data. As a result of extensive data 
comparisons, the model listed below for eddy viscosity is believed to adequately incorporate the 
effect of compressibility and provide a universal description over a range of Mach numbers: 


€ 


( 


Ue, 

P 

^ U r K Y, 


Y > Yj, 

Y < Y x , 


(15) 


where Y x = p e 8* U e K//c. 


This is the modification of the Cebeci-Smith model while, for the 
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Baldwin-Lomax model, the first of equations (15) is replaced by 


e 


P max P w 

P* 


C P K Y„ 


( 16 ) 


In these equations, subscripts w and e denote quantities evaluated at the wall and in the 
mainstream respectively. A similar model has also been obtained for heat transfer and this, as 
well as a detailed description of a portion of the work, is described in the accompanying paper. 

An an example of the type of profile produced in the present analysis, the predicted 
velocity profile using the model in equations (15) and (16) is shown in Figure 1; the data was 
taken at a Mach number of 4.5 for an adiabatic wall. Note that this is not a curve fit and, 
beside the wall and mainstream temperature and velocity, the only input is <5*, the 
incompressible displacement thickness at that streamwise location. In Figure 2 a 
corresponding plot is shown, where the conventional “incompressible” Baldwin-Lomax model 
(given in equation (7)) has been used. It is clear that an “incompressible” model is 
inadequate; in additon, it emerges that the skin friction is considerably overpredicted by the 
“incompressible” model. Further comparisons over a Mach number range up to 10 are given 
in the attached paper. 

To date, some important aspects resolved by the research program are as follows: 

(1) The use of the Howarth-Dorodnitsyn variable 

Y = [^ pdy 
J 0 


is critical. Both the velocity and total enthalpy profiles exhibit a logarithmic behavior in the 
overlap zone over a wide range of Mach numbers. The profile behavior appears to be universal 
when expressed in terms of Y. 

(2) Standard turbulence models which have been tuned for low Mach number flows 
should not be used in an unaltered fashion as the Mach number increases. The modifications 
introduced in the present study involve factors of the density and appear to be universal with 
increasing Mach number. 
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(3) A new model has been developed for heat transfer in a high-speed compressible 
flow. The model was developed independent of a Reynolds analogy argument or a turbulent 
Prandtl number concept. The correspondesnce with existing measured total enthalpy data is 
very encouraging. 

(4) One byproduct of the present research is a set of self-similar profiles for total 
enthalpy and velocity which can be used to initiate a Navier-Stokes calculation, simply by 
estimating a displacement thickness distribution along each wall. 

(5) A set of embedded functions for total enthalpy and velocity has been developed for 
supersonic flows. 

3. Overview of UTRC Computational Work 

The following discussion summarizes the work being carried out at UTRC by Greg D. 
Power and Seyhan Ersoy. 

The analytical wall layer model of Walker 1 has been incorporated into the NASCRIN 
full Navier-Stokes (NS) code developed originally by Kumar 2 . This model provides an 
analytical description of the velocity and energy profiles in the inner layer of a turbulent 
boundary layer. Therefore, the numerical solution technique for the bulk of the flow need not 
be applied in the region where the wall layer model is applicable. In order to accurately 
determine skin friction and heat transfer, standard solution techniques require up to half of the 
grid nodes be placed within the inner layer. Since the wall layer model replaces the numerical 
solution in the inner layer, computer storage is reduced and, indirectly, the CPU time required 
to reach a steady-state solution is also reduced. As a preliminary step, this concept has been 
incorporated into the ABLE 3 boundary-layer code by Walker et al 4 . In the following 
discussion, a review of the wall layer concept in the ABLE code, and several aspects of the 
implementation of this concept in NASCRIN, are presented. 

Algorithm 

The ABLE code is a direct boundary-layer solution procedure. The governing 
equations were originaly written using a turbulent flow extension of the Levy-Less 
transformation. In this transformation, the normal coordinate is scaled by the zero-pressure 
gradient boundary-layer growth resulting in a transformed coordinate, r /. In the wall-layer 
model version of ABLE, the normal coordinate is scaled by a parameter A. This parameter is 
adjusted as a function of the solution at each marching step in order to maintain y + = 60 at 
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the first grid point from the wall, i.e. the grid is adapted to the edge of the inner layer. The 
points which would normally be needed within the inner layer are deleted, thus reducing the 
number of points across the entire boundary laeyr by approximately half. 


The critical information used to determine the solution at the second grid point (^ 2 ) is 
that the outer defect law and the inner layer solution match in the log-law region 
(50 < y + < 100). In this region, the inner and outer velocity profiles take the form: 

Inner Solution: y + —» oo 


U + = \ log y + + c,-, (18) 

Outer Solution : t? — + co 

F = 1 + log i) + Co), (19) 


where 


+ _ £ 


u = 


F = jJL and 


u e > 


Also the independent variables are given by 


+ _ fhu UtY 

y — & w 



where the definition of the normal coordinate will be discussed later. 

The unknowns in these equations are u*, C 0 , and F 2 (the solution at tj 2 ). A third 
condition is required in order to close the system. This condition is that the solution at t] 3 
(F 3 ) must also obey the log law given by equation (19). The solution at t] 3 is a result of 
integrating the boundary-layer equations across the layer using, in this case, the Cebeci-Smith 5 
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turbulence model. The boundary-layer solution and the solution at tj 2 are determined 
iteratively at each streamwise station in two steps. 

1. For an assumed u* and F 2 , equation (19) is solved at tj 2 for C 0 . Then, the velocities 
from equation (1) and (2) are matched to determine a new u*, i.e. 


Iu = K lo S 


PujMt a 


+ C t * - Co, 


2. Equation (19) is evaluated at and r} 3 , then these equations are subtracted to get 



The boundary-layer equations are solved with this equation as the lower boundary 
condition. 


The above steps are repeated until the solution at the current marching station converges. In 
summary, the value of u* is determined from matching the solution of to the known inner 
layer solution and the flow field is determined using a gradient boundary condition which is 
based on the log law behavior of the outer layer veloctiy profile. Note that a similar iteration 
scheme is used for enthalpy. 

It is not practical to adjust the NS grid in order to adapt to the inner layer thickness. 
Since the only requirement for the match point is that this point lies within the log-law region, 
a fixed grid can be used where the solution below y + w 60 is evaluated analytically and the 
solution above the match point is computed using the numerical algorithm with a boundary 
condition set at the match point. In this matter, the grid location of the match point can 
change as a function of time and space. 

NASCRIN is an explicit time-marching procedure. The boundary conditions are set at 
the boundaries for each time based on the solution in the flow field at the previous time. In 
the original formulation, the wall layer boundary condition is implemented at the match point 
in a two-step procedure: 


10 



1. The friction velocity, U T , is evaluated from the solution at the match point and the 
first point above the match based on the log law behavior of the outer profile: 



2. The solution at the match point is updated by matching the analytical inner layer 
solution (u + ) and the outer solution from NASCRIN: 


llm 


u T u + (ym). 


If the solution is not converged, the friction velocty computed in step 1 could force the 
solution at the match point (step 2) to change substantially. Since the explicit time-marching 
algorithm does not react immediately to changes in boundary conditons, the next time step 1 
is implemented and there could be a large change in u r . Repeating these steps as the solution 
progresses eventually results in the solution diverging. Two methods to stabilize the iterative 
scheme are to under-relax the u r calculation or to implement the wall layer model every few 
time steps rather than every time step. Both of these techniques slowed, but did not 
eliminate, the divergence. A third technique was attempted which is related to experimental 
methods for determining the skin friction from velocty profiles. Instead of using the solutions 
at m and m+i in step 1, the gradient of the log plot is determined over the entire range for 
which the log law is applicable. Therefore, the solution at the match point does not have as 
much of an influence on the evaluation of u T . This technique also slowed, but did not stop, 
the divergence. 

Note that the steps outlined above for NASCRIN are implemented in the reverse order 
to the ABLE implementation, i.e. the solution gradient is used to determine u T and the inner- 
layer analytical solution is imposed as a boundary condition. Therefore, a new algorithm was 
constructed for the NS code which more closely resembles the ABLE algorithm: 

1. u r is determined from the solution at one point above the match point by matching 
the analytical inner layer solution (u + ) and the outer solution from NASCRIN: 
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m + 1 

+ (ym + l) 


2 . 


The solution at the match point is determined through an explicit implementation of 
the gradient boundary condition based on the log law behavior of the outer solution: 


Um = U m+1 + ^ log 


Y 


m+1 


There are two advantages to this approach. First, u r does not vary significantly from 
iteration to iteration since the solution at m+i slowly adjusts to the new boundary condition. 
Second, the gradient boundary condition is somewhat weaker so that the boundary condition 
also does not vary rapidly. Preliminary results using this technique have been encouraging. 

Compatibility 

In order for the above algorithms to succeed, the outer solution provided by the 
boundary-layer or the Navier-Stokes procedure must be compatible with the inner-layer 
anlaytical solution. The value of the friction velocity is determined from the Navier-Stokes 
solution based on the analytical inner layer solution. The friction veloctiy is then used to set 
the boundary condition at the match point. If the velocity profile in the region where 
NASCRIN is applied does not behave identically to the inner layer solution, the value of u r 
will be in error and thus the entire solution will be incorrect. The turbulence model is the 
primary influence on the velocity profile. 

Most turbulence models, including the Baldwin-Lomax model, were developed for 
incompressible flows. At supersonic conditions, the density can vary greatly across the 
boundary layer and, thus, the incompressible assumption is invalid. By correlating with 
numerous experimental data sets, Walker 6 has found that the normal coordinate (Y) should be 
defined as 
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Using this definition, the logarithmic regions of u + vs. y + plots collapse to a single 
curve independent of Mach number and wall temperature. Traditionally, Y has been set to 
the physical distance (y) from the wall. Using this incompressible formulation can lead to 
large errors in eddy viscosity. For example, the density at the wall of an insulated flat plate 
immersed in a Mach 2.5 flow is nearly one-half of the free-stream density. Using the 
compressible formula, the eddy viscosity in the inner layer can be shown to be 

Hi = *| K 2 Y 2 gj for y < y*. 


The Baldwin-Lomax 7 outer layer eddy viscosity takes the form 


f*0 


PePw 

P 


KCcpymaxFmar 


for y > y c , 


where y c is the cross-over at which /i, < /i 0 . 

Prog res s 

The algorithm described above for incorporating the wall-layer model into the 
NASCRIN code requires that boundary conditons be set at a point within the flow field and 
that the solution below this point not be updated. First the wall-layer boundary conditions for 
the streamwise velocity and enthalpy are set at the match point in a new subroutine, 
WALFUN. In this subroutine, the analytical solutions for these variables are also set from the 
match point to the wall. Two other boundary conditions are required. The normal velocity at 
and below the match is determined by integrating the steady-state continuity equation. Since 
the log-law region is well within the boundary layer, the normal pressure gradient must be 
zero, and, thus, the pressure within the inner layer is set using this fact. 

In subroutine VISCOS, where the explicit integration takes place, the time step (At) is 
set to zero at and below the match point. By doing so, the solution in the inner layer is not 
updated in the integration scheme. Note that the match point need not be a constant grid line 
as a function of time or streamwise distances. 
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The logic described above has been implemented in NASCRIN. In order to check the 
algorithm, first a converged solution of NASCRIN with integration to the wall was obtained, 
A second run was made by simulating the wall-layer implementation, i.e. the converged 
solution for streamwise velocity and enthalpy at the match point was held fixed and the other 
variables were set as described above. The solution converged to essentially the same solution 
as when integrating to the wall. 

The compressibility corrections for the Baldwin- Lomax turbulence model have also 
been incorporated. The code has been applied using the original and modified models to a test 
case corresponding to the experimental configuration studied by Coles 8 . The configuration 
consists of an insulated flat plate immersed in a Mach 2.5 flow at a Reynolds number of nearly 
3 x 10 6 /meter. Skin friction and a limited number of velocity profiles are available for 
comparison. The results computed using the original turbulence model overpredict the skin 
friction data while the results computed using the modified model underpredict the data as 
shown in Figure 3. This study indicates that there is a significant impact of the 
compressibility corrections on the solution. 

Implementation of the wall-layer model is currently being studied. Preliminary results 
indicate that the second formulation described in the algorithm section of this memo is stable 
and converges. The time step required for stability of the explicit integration scheme was 
increased due to the wall-layer implementation by a factor of almost 2.5 compared to the full 
NASCRIN solution. However, the results using the wall-layer model do not match the data or 
the NASCRIN solution when integrated to the wall. A comparison of u r computed with and 
without the wall-layer model is shown on Figure 4. The primary discrepancy is that the slope 
of the u r (skin friction) curve does not match the expected slope. The compatibility of the 
inner-layer analytical solution and the outer-layer NASRIN solution may be the source of the 
error. In fact, the discrepancy in this solution is most likely related to the discrepancy found 
with the turbulence model above. 

At this time, the algorithm is stable and leads to converged solutions. The 
compatibility of the inner-layer solution with the turbulence model is being studied. 
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Figure 1. Comparison of velocity data with the theoretical profile using the new compressible 
Baldwin- Lomax model. 





Y + 


Figure 2. Same data as in Figure 1, but using the conventional Baldwin-Lomax Model. 
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Figure 3. Effect of turbulence model compressibility on NASCRIN solution. 
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